This notebook is part of the HYDRA-EO project (ESA, CROP
MULTIPLE STRESSORS). Here we demonstrate how to use the
project’s in-house R packages
(ToolsRTM) to simulate canopy reflectance
(BRF), and vegetation indices (NDVI).
The aim is to illustrate how radiative transfer
models can be used to
generate synthetic datasets under different stress
scenarios (biotic/abiotic),
which are then used to validate remote sensing
algorithms and analyze the
sensitivity of spectral indices to changes in vegetation traits.
References:
- HYDRA-EO is funded by the European Space Agency (ESA) under the EXPRO+ Tender “Crop Multiple Stressors, Pests and Diseases”. Action 1-12684
- Camino et al., ToolsRTM package (GitLab: ToolsRTM)
- Camino et al., SCOPEinR package (GitLab: SCOPEinR)
To begin the exercise, we ensure that all necessary R packages are
installed and loaded.
These libraries support a range of functionalities, including:
ggplot2,
plotly, leaflet, leafem,
htmlwidgets)terra, sf, stars,
gdalcubes)dplyr,
tidyr, DT, magrittr)rstac,
fs)MASS,
fields,corrplot)viridisLite, RColorBrewer,
leaflet.extras, reshape2)The script below checks for any missing packages and installs them if needed, ensuring a smooth experience when running the rest of the notebook.
Install the ToolsRTM Package
The ToolsRTM package is an R toolkit developed for
canopy radiative transfer modeling. It integrates forward simulations
from well-known models such as PROSPECT (leaf),
SAIL (canopy architecture), and their coupled form
PROSAIL.
With these RTM models you can:
This package is hosted on GitLab. It is
not on CRAN, so we install it directly from GitLab
using the remotes package. For exploring more options, you
can open the documentation in your browser:
browseURL("../docs/ToolsRTM.html")
After installation, you can simply load the package with
library(ToolsRTM)in future sessions.
if (!requireNamespace("ToolsRTM", quietly = TRUE)) {
if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
# Public GitLab repo (no token needed). Set upgrade="never" for reproducibility.
remotes::install_gitlab("caminoccg/toolsrtm", upgrade = "never")
}
library(ToolsRTM)
##
## Attaching package: 'ToolsRTM'
## The following object is masked from 'package:utils':
##
## fix
Before running systematic simulations in R, we will use the interactive Shiny app to explore how changes in leaf biochemistry and canopy structure affect reflectance spectra. This step provides intuition for interpreting results later.
Launch the online simulator (interactive):
The ToolsRTM package includes the helper function
getSim_fromLUT(), which allows us to vary a single input
trait across a chosen range, keep all others fixed, and directly
visualize the resulting spectra (using
method = "ggplot").
This is a simple way to see how sensitive canopy reflectance is to
changes in one parameter at a time.
Below, we test one examples:
# Cab from 10 to 80 µg/cm², step 10 (Interval)
cab_plot <- ToolsRTM::getSim_fromLUT(
trait = "Cab",
nmin = 10, nmax = 80, Interval = 10,
model = "PROSAIL",
method = "ggplot"
)
In addition to Cab, the inputs object
in ToolsRTM defines many other traits that can be
varied in getSim_fromLUT() or when building LUTs. These
represent both leaf-level biochemistry/structure and
canopy/scene parameters:
Leaf traits
Car — carotenoids (µg/cm²)Anth — anthocyanins (µg/cm² or relative units)Cbrown — brown pigments (0–1, senescence/stress)LMA — leaf mass per area (g/cm², replaced by
Prot + CBC in PROSPECT-PRO)EWT — equivalent water thickness, same as
Cw in many PROSAIL formulations (g/cm²)Prot — protein fraction (g/cm², part of LMA in
PROSPECT-PRO)CBC — carbon-based constituents (cellulose + lignin;
g/cm²)N — leaf structure parameter (dimensionless; controls
scattering)Canopy traits
LIDFa, LIDFb — leaf angle distribution
parameters (Campbell distribution)hspot — hotspot parameter (0–1, controls BRDF hotspot
intensity)LAI, — leaf area index (m²/m²)Scene geometry
tts — solar zenith angle (°)tto — sensor/view zenith angle (°)psi — relative azimuth angle (°)In this section you will:
Build the LUT (100 samples)
We start from the default PROSAIL input parameter ranges provided in the ToolsRTM package. A LUT is then generated by random sampling across these ranges, in this case with 100 samples. This LUT will be the basis for running the simulations in forward mode.
#Get default PROSAIL input list
inputs <- ToolsRTM::inputsPROSAIL
# (Everything else stays as in defaults; you can tweak similarly if needed)
# 2.2 Generate a LUT with 100 samples
set.seed(1234)
LUT <- as.data.frame(ToolsRTM::getLUT(inputs = inputs, nLUT = 100, setseed = 1234))
# Summarise min–max per parameter
lut_ranges <- data.frame(
Parameter = names(LUT),
Min = round(sapply(LUT, min, na.rm = TRUE),3),
Max = round(sapply(LUT, max, na.rm = TRUE),3)
)
# Display in an interactive datatable
DT::datatable(
lut_ranges,
caption = "Ranges of sampled input parameters in the LUT",
options = list(pageLength = 10)
)
In PROSAIL, the soil spectrum is a critical boundary condition, especially when the canopy is sparse (low LAI). Soil reflectance can strongly influence the visible and near-infrared regions of canopy reflectance.
The ToolsRTM package provides reference soil spectra:
dry_soil → brighter across VIS–NIR, typical of bare or
dry soils.wet_soil → darker and flatter, since water reduces soil
reflectance.A mixed soil spectrum is created as a weighted combination of the
two, controlled by psoil:
rsoil.dry <- ToolsRTM::dataSpec_PDB$dry_soil
rsoil.wet <- ToolsRTM::dataSpec_PDB$wet_soil
psoil <- 0.5
rsoil.default<- c(psoil*rsoil.dry+(1-psoil)*rsoil.wet)
# Data frame for ggplot
soil_default_df <- data.frame(
Wavelength = 400:2500,
Reflectance = rsoil.default
)
# Plot with ggplot
p.soil <- ggplot(soil_default_df, aes(x = Wavelength, y = Reflectance)) +
geom_line(color = "blue", size = 0.8) +
labs(
title = paste0("Default Mixed Soil Spectrum (psoil = ", psoil, ")"),
x = "Wavelength (nm)",
y = "Reflectance"
) +
theme_bw()
# Convert to interactive plotly
ggplotly(p.soil, tooltip = c("x","y")) %>%
config(displaylogo = FALSE)
For each LUT entry, we run the PROSAIL fourSAIL canopy model using PROSPECT-PRO at the leaf level. We assume a mixed soil spectrum (50% dry, 50% wet). The result is a BRF spectrum for each of the 100 LUT samples.
wl_grid <- 400:2500 # wavelength grid in dataSpec_PDB
#Simulate PROSAIL (foursail) for each LUT row
sims <- lapply(seq_len(nrow(LUT)), function(i) {
out <- suppressMessages(ToolsRTM::foursail(
inputLUT = LUT[i, ],
rsoil = rsoil.default,
LeafModel = "PROSPECT-PRO",
spectrum.all = TRUE
))
# out has rdot, rsot, rddt, rsdt; combine to BRF using solar zenith (tts)
brf <- ToolsRTM::Compute_BRF(
rdot = out$rdot,
rsot = out$rsot,
tts = LUT[i, "tts"],
data.light = ToolsRTM::dataSpec_PDB
)
# Return tidy tibble for this run
tibble(run = i, wl = wl_grid, rho = as.numeric(brf))
})
Once all simulations are complete, we aggregate the results to compute a mean spectrum ± standard deviation. This shows the overall spectral variability caused by the input parameter ranges.
# 1) Combine all runs
brf_runs <- dplyr::bind_rows(sims)
# 2) Summarise mean ± sd by wavelength
spec_stats <- brf_runs %>%
group_by(wl) %>%
summarise(mean_rho = mean(rho, na.rm = TRUE),
sd_rho = sd(rho, na.rm = TRUE),
.groups = "drop")
The plot below shows all 100 PROSAIL simulations
(gray lines), the mean spectrum (blue line), and the
variability envelope (shaded area = ± 1 standard
deviation). This allows us to see which spectral regions are most
variable across the sampled LUT.
For example:
- In the visible region (400–700 nm), variability is
driven mainly by pigment content (Cab, Car, Anth, Cbrown).
- In the NIR plateau (700–1300 nm), variability
reflects structural and canopy parameters such as LAI and leaf angle
distribution.
- In the SWIR (1400–2500 nm), differences in water
(EWT) and dry matter (Prot, CBC) dominate.
Explore the interactive plot by hovering to read wavelength–reflectance values and toggling the legend items.
p <- ggplot() +
geom_line(data = brf_runs, aes(x = wl, y = rho, group = run), alpha = 0.12) +
geom_ribbon(data = spec_stats,
aes(x = wl, ymin = mean_rho - sd_rho, ymax = mean_rho + sd_rho),
inherit.aes = FALSE, alpha = 0.2, fill = "blue") +
geom_line(data = spec_stats, aes(x = wl, y = mean_rho),
inherit.aes = FALSE, color = "blue") +
labs(title = "PROSAIL (custom LUT): mean ± 1 sd (BRDF)",
x = "Wavelength (nm)", y = "Reflectance (BRDF)") +
theme_bw()
# Convert to interactive plotly
ggplotly(p) %>%
config(displaylogo = FALSE)